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ANALYSIS OF HEAT TRANSFER IN A BOUNDED TWO-DIMENSIONAL POROUS 
REGION AND APPLICATION TO A POROUS SQUARE DUCT 
by Marvin E. Goldstein and Robert Siegel 
Lewis Research Center 

SUMMARY 

A finite two-dimensional permeable porous region is considered bounded by four 
surfaces. Two of the opposing surfaces have no fluid flow or heat transfer across them 
and are consequently along the direction of the fluid and heat flows. The other two sur- 
faces are each in contact with a reservoir at a different constant pressure, and this pro- 
duces a flow toward the low pressure reservoir. The flow is assumed to be governed by 
Darcy's law, and hence, the pressure acts as a velocity potential. As a result of the 
boundary conditions the porous region occupies a rectangular region in a complex poten- 
tial plane. The energy equation for the temperature distribution in the medium is shown 
to be transformed into a separable equation when its independent variables are changed 
to the coordinates of the potential plane. The general solution is obtained in the potential 
plane for the cases where the side of the porous material at which the flow exits is either 
maintained at an arbitrary temperature or has an arbitrary imposed heat flux. Confor- 
mal mapping is then used to relate the solution to the physical geometry. The method is 
illustrated by obtaining the heat -transfer characteristics of a square duct with porous 
walls of finite thickness. 


INTRODUCTION 

A method for extending the use of a metallic structural material to higher temper- 
ature applications is to provide transpiration cooling. The metal is made in a porous 
form , and coolant is forced through it from a reservoir toward the boundary exposed to 
the high temperature source. Some possible applications are for cooling turbine blades, 
rocket nozzles, arc electrodes, and portions of surfaces during either high speed flight 
or reentry into the Earth's atmosphere. 

The energy equation governing the temperature distribution in the permeable porous 



material contains the velocity distribution. The velocity is generally a complicated 
function of position, which makes it difficult to solve the energy equation analytically. 

The solutions in the literature have consequently been limited to one -dimensional situa- 
tions. However, in many applications, the geometries are two- or three-dimensional. 
Solutions can be obtained numerically, but it would be desirable to have some analytical 
results to check numerical and approximate procedures. In addition, formulas derived 
from analytical solutions can be used for calculating results at the porous boundaries 
without obtaining detailed calculations for the interior of the region as in a numerical 
solution. This is convenient for coupling the heat transfer in the porous medium into the 
overall heat transfer system such as coupling it with an external thermal boundary layer. 
Analytical results are also useful in evaluating the importance of the governing param- 
eters and in helping to determine conditions when two- or three-dimensional effects be- 
come significant so that a locally one -dimensional solution will not apply. 

For these reasons the present authors devised in reference 1 an analytical method 
for obtaining the heat -transfer behavior in a two-dimensional porous medium. The the- 
ory was developed for a porous wall of arbitrarily varying thickness but of infinite length. 
The purpose of the present report is to extend the theory to finite two-dimensional re- 
gions. The example used to illustrate the general solution that is obtained is the heat 
transfer in a square duct with porous walls of a finite thickness. 

The following is a brief outline of the analytical method derived in reference 1 that 
will be further developed and applied herein. As stated in reference 1, it will be assumed 
that the fluid and solid matrix are locally in good thermal contact so that the local fluid 
and solid matrix temperatures are equal. As a consequence, a single energy equation 
can be written that includes the heat transport by conduction in the matrix as well as by 
coolant convection. 

The velocity in the convective term of the energy equation is a function of position 
and is proportional to the local pressure gradient for the slow viscous flow often en- 
countered in porous media. The pressure can therefore be regarded as the velocity po- 
tential. The equations governing the problem are transformed into a potential plane. 

Since two boundaries of the porous medium are at constant pressure, these become par- 
allel lines in the potential plane. The other two boundaries have no flow across them and, 
hence, are perpendicular to the constant potential boundaries. Thus, in the potential 
plane the porous region becomes a rectangle. 

Since the geometry in the potential plane is simply a rectangle regardless of the 
physical shape of the porous material, a solution in this region of the potential plane will 
apply to all geometries. An important factor is that, when the energy equation is trans- 
formed into the potential plane, it becomes a separable equation, and hence, the general 
solution can be obtained for the rectangular region. This solution can then be trans- 
formed into specific physical geometries by conformal mapping. If the geometry is too 
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complicated to map analytically, the mapping can be performed numerically; one tech- 
nique for this is given in reference 2. 
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SYMBOLS 


half -width of outside of square duct divided by wall thickness 
half -width of outside of square duct 
coefficient defined following eq. (48) 
specific heat of fluid 
incomplete elliptic integral of first kind 
function in specified heat flux distribution 
function in specified temperature distribution 
reference length in porous material 

complete elliptic integral of first kind; K'(k) = K^l - k^ 
thermal conductivity of porous region 
quantities defined in eq. (A8) 

dimensionless coordinate along boundary S, i f h 

b X 

coordinate along boundary s measured from left side 
heat flux parameter (qg - q^)/qj 
temperature parameter (tg - t^)/(tj - tj 
unit outward normal vector 
pressure 
heat flux 


bounding surfaces of porous region in dimensionless coordinate system 
bounding surfaces of porous region 

dimensionless temperature defined in eqs. (11) and (12) 
temperature 

a h r _ 

dimensionless velocity, it i u 

K Po - P S 


velocity 
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u intermediate mapping plane, a + i/3 
V dimensionless velocity component in Y-direction 

v velocity component in y-direction 

W complex potential, ip + icp 

X,Y dimensionless coordinates, x/h r and y/h r 
x,y rectangular coordinates 

Z dimensionless physical plane, X + iY 

a real part of u -plane 

/3 separation constant in solution for ®; imaginary part of u-plane 

y n quantity defined in eq. (48) 

A nondimensionalizing quantity defined by eq. (12) 

7] imaginary part of co 

® variable defined in eq. (38) 

k permeability of porous material 

PC D K K ' Pfl> 

X parameter, — £_ — _ 

2 k p 

m 

M fluid viscosity 

£ real part of co 

p fluid density 

<f> function of w in solution for ® 
cp potential (p n - p)/ (p^ - p_) , imaginary part of W 
function of ip in solution for © 
ip real part of W 

n quantity defined in eq. (A15) 

co intermediate mapping plane, £ + irj 

V dimensionless gradient defined in eq. (18) 

Subscripts: 

l at left boundary 

o at lower boundary 
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r at right boundary 

s at upper boundary 

1, 2 values on upper boundary at left and right sides 

oo coolant reservoir 

GENERAL ANALYSIS OF TWO-DIMENSIONAL POROUS COOLED REGION 

Governing Equations 

Consider the two-dimensional porous region with effective thermal conductivity k m 
(based on the entire cross sectional area) and permeability k shown in figure 1(a). The 
lower surface of the region whose unit outward -drawn normal is n Q is denoted by s Q 
and the upper surface of the region whose unit outward -drawn normal is n g is denoted 
by s. The left and right surfaces of the region whose unit outward -drawn normals are 




(b> Dimensionless physical plane. 

Figure 1 - Two-dimensional porous region with no fluid or heat flow through two opposite 
sides. 
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and n r , respectively, are denoted by and s r , respectively. There is a fluid with 
constant density p, constant heat capacity C p , and constant viscosity ju i which is flow- 
ing through the region. Assume that the thermal conductivity of the fluid is very small 
compared with k m and that the pore size is so small that Darcy's law holds. Let u 
denote the Darcy velocity (local volume flow divided by entire cross section rather than 
by open area) of the fluid. We suppose that no changes occur in the direction perpen- 
dicular to the x,y-plane so the situation is two-dimensional. 

Below the porous region (see fig. 1(a)) there is a reservoir which is maintained at 
constant pressure and temperature p Q and t^, respectively. There is no heat or mass 
flow across the surfaces s^ and s r - The pressure of the fluid above the wall is con- 
stant and equal to p g . We suppose that p Q > p g . Then the fluid flows from the reser- 
voir through the porous region and out through the top surface. Since p rt and p are 

o s 

both constant, the fluid velocity at both the upper and lower region surfaces will be in a 
direction perpendicular to these surfaces. 

If the thermal communication between the fluid and the matrix is sufficiently good, 
the local fluid temperature will be approximately equal to the local solid matrix temper- 
ature. We denote this common temperature by t. When these assumptions are made, 
the heat and mass flow within the porous material are governed by the equations given in 
reference 1. Thus, 


k m V 2 t - PCpU • Vt = 0 (1) 

u = - — Vp (2) 

V ■ u = 0 (3) 

Boundary Conditions 

As the fluid in the reservoir approaches the porous region, the fluid temperature 
will rise from the reservoir temperature t^ to the region surface temperature t 
which is an unknown variable along s Q . Since the thermal conductivity of the fluid is 
assumed to be much less than the thermal conductivity of the matrix material, the thick- 
ness of the fluid layer over which this temperature rise takes place is very small com- 
pared with the porous region thickness provided the flow is not very small. We can, 
therefore, assume that the liquid layer is locally one-dimensional. Since the velocity is 
perpendicular to s q , there is no flow along this thermal layer. Hence, applying an 
energy balance to the thermal layer shows that 
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(4) 


k m% ' Vt = pC p (t " too) % ' H 


P = P r 


for (x,y)es Q 


Since there is no flow of heat or mass across the side surfaces s^ and s r , the normal 
derivative of both the pressure and the temperature must vanish at these surfaces. 
Hence, 


n l - Vt = 0 > 

► 

n z • Vp = 0 > 


for (x,y)es z 


n r • Vt = 0 ^ 
h r • Vp = 0 ^ 


for (x,y)es r 


(5) 


( 6 ) 


We shall consider two different types of thermal boundary conditions for the upper sur- 
face s. First, we shall suppose that its surface temperature is specified so that it 
varies along the surface from a temperature t^ at the left end to a temperature t 2 at 
the right end. Let l be the distance measured along the upper surface from its left 
end. Then for this case the boundary conditions on the upper surface of the region are 


‘ = l s = ‘i + <‘ 2 - l i> H l r 

n ry 


> for (x,y)es 


(7a) 


P = P C 


J 


where H is a given function equal to zero at (Z /h y ) = 0 and equal to 1 when (£ g /h r ) is 
a maximum. The constant h r is an appropriate reference length. 

We shall also consider the case where the heat flux into the wall is specified along 
the upper surface and it varies from a value of q^ at the left end to a value of qg at 
the right end. Under these conditions the boundary conditions on the upper surface be- 
come 
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> for (x,y)es 


(7b) 


k m”s * Vt = q g = qj + (q 2 - q^ 

P-P s 

where G is a given function which is equal to zero at l / h equal to zero, and equal to 

S X 

1 when l / h is a maximum, 
s r 

Equations (1) to (3) together with either boundary conditions (4) to (6) and (7a) or al- 
ternatively the boundary conditions (4) to (6) and (7b) completely determine the solution 
to the heat transfer and flow problems within the porous wall. 



Dimensionless Equations 


It is now convenient to introduce the following dimensionless quantities: 


' . PC p "ftp - p s> 
2k m " 


N = 


t 2 ~ t l 

t 1 - t 

1 °° 


M = 


*2 -*1 


( 8 ) 

( 9 ) 


( 10 ) 


where 




X = — 


Y = y. 

h_ 


L s = — 
h r 


(p = 


Pq -P 

Po-P S 




u = £ _ — - — u 

K (p 0 - p s ) 


t - 1 


T = 


A 


J 


(ID 


(tj - t^) if boundary condition (7a) applies 


AM 


^l h r . 


if boundary condition (7b) applies 


m 


(12) 


Upon substituting these definitions into equations (1) to (3) and boundary conditions 
(4) to (7b), we obtain 


v T - 2A.U • VT = 0 
U = Vcp 
V • U = 0 


(13) 
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(14) 


where 


n • VT = 2An 0 • UT~^ 


cp = 0 


■ VT = 0 


• Vcp = 0 


n r ■ VT = 0 


n r • V (p = 0 


for (X, Y)eS 


for (X, Y)eS, 


for (X, Y)eS 


(15) 


(16) 


T = 1 + NH(L g ) 


for (X, Y)eS 


(17a) 


or 


(p = 1 


n • VT = 1 + MG(L.) 

S o 


cp= 1 


for (X, Y)eS 


(17b) 


^ r) r) 

V 

ax aY 


(18) 


and the porous region in dimensionless coordinates is shown in figure 1(b). 

The second equation (13) can be used to eliminate U in the other two equations (13). 
Thus, 


V 2 T - 2A Vcp • VT = 0 


(19) 


and 
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( 20 ) 


~2 » 
V </> = 0 


Since cp is constant on both S Q and S, it is clear that 


n_ = 


0 M 


for (X,Y)eS 


(21a) 


and 


n = -*2- 


V qj 


for (X, Y)eS 


(21b) 


Using these results together with the second equation (13) in the boundary conditions 
(14) and (17b), we obtain 


and 


Vcp • VT = 2XT |V<p|‘ 
& = 0 


> for (X,Y)eS 


(22a) 


V0 ■ VT = |Voj| [l + MG(L S )] 

(P = 1 


> for (X, Y)eS 


(22b) 


Porous Region in Potential Plane 

Since equation (20) shows that tp satisfies Laplace’s equation, there must exist a 
harmonic function tp and an analytic function W of the complex variable 


Z = X + iY 


such that 


W = ip + i<p 


(23) 


(24) 
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Physically the change in ip between any two points is proportional to the volume flow of 
liquid crossing any curve joining those two points. Hence, since there is no flow across 
either surface or S r , the function ip must be constant on both of these surfaces. 
Since W is determined only to within an arbitrary constant, we can always arrange 
matters so that 


4> = 0 for (X, Y)eSj 

We shall denote the constant value of ip on S r by ip r - Thus, 


9mH = 0 

for 

ZeS ( 

= 1 

for 

ZeS 

= 0 

for 

ZeS z 

Re W = ip r 

for 

ZeS r 


This shows that the mapping 


Z - W 


(25) 


transforms the region occupied by the porous material in the dimensionless physical 
plane (Z -plane) shown in figure 1(b) into the rectangular region in the W-plane shown 

v 



Figure 2. - Porous region in complex potential plane W = ijj + \<p. 
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in figure 2. The boundaries S Q and S in the physical plane are transformed into the 
lines cp = 0 and 0 = 1 , respectively, in the W-plane, and the boundaries S 7 and S 
in the physical plane are transformed into the lines ip = 0 and ip = ip^,, respectively, in 
the W-plane. The mapping Z — W can be found by using the techniques of conformal 
mapping once the shape of the porous wall in the physical plane is specified. This gives 
the solution to the boundary value problem for the mass flow. (An example for a spec- 
ific wall geometry will be worked out subsequently.) 

Since ip is constant on both and S r , it is clear that 


n 7 



for (X,Y)eS z 


for (X,Y)eS 

N | 


Using these results together with equation (25) in the boundary conditions (15) and 
(16), we obtain 


and 


Vip • VT = 0 

ip = 0 


for (X,Y)eS z 


(26) 


Vi p ■ VT = (T 

>■ 

* = r , 


for (X,Y)eS r 


(27) 


Transformation of Boundary Value Problem Into Potential Plane 

The boundary value problem for the temperature in the physical plane can be solved 
if we use the Boussinesq transform (ref. 3) to transform it into a boundary value prob- 
lem in the W-plane. Thus, the independent variables X and Y in equation (19) and 
its boundary conditions will be changed to the variables 0 and ip, and the resulting 
boundary value problem will be solved in the rectangular region in the W-plane shown 
in figure 2. Once this solution which gives T as a function of <p and ip has been found, 
the mapping W — Z (which is completely determined once the geometry of the region in 
the physical plane is specified) can be inverted to give Z as a function of 0 and ip. 
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Thus, T will be known parametrically (in terms of the parameters <p and ip) as a func- 
tion of X and Y. This will complete the solution to the problem. 

Apropos of these remarks, recall that (ref. 4) 


V 




+ 




and 


Notice that 


r^j 

Vcp 


dW 

dZ 


(28) 


(29) 


dW 

_ dW 

/dZ dw\* = 

dW 1 2 /d 

\ * 

dZ 

dZ 

\dW dZ / 

dZl \d 

W/ 


Hence, upon taking real and imaginary parts, we find 


dcp _ 

dW 

2 ax 

0X 

dZ 

dcp 

3 cp „ 

dW 

2 3Y 

9Y _ 

dZ 

dcp 


and 


d ip_ = 

dW 

2 ax 

ax 

dZ 

dip 

dip_ _ 

dW 

2 3Y 

dY 

dZ 

dip 


This shows that 
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2 


(30) 


Vcp • VT = ^iX + ST Ml 

dW 

2 

9T 

dW 

\3X dcp 3Y 9<p/ 

dZ 


dcp 

dZ 


and 


/V 

Vip ■ 




(31) 


Finally, notice that since cp is constant on S, the distance L g along S is a function 
of ip only. (This functional relation is known once the mapping W — Z which solves 
the flow problem is known. ) Thus, 


L S (X,Y) = L s (ip) for (X,Y)eS 


(32) 


Now, using equations (28) to (32) in equation (19) and the boundary conditions (22a), 
(26), (27), (17a), and (22b) yields 


or 


+ 2x11 = 0 

dip 2 dcp 2 d(P 


— - 2AT = 0 for cp = 0 

dcp 


— = 0 for ip = 0 

dip 


— = 0 for ip - ip 

r 


T = 1 + NH(L g (^)) for cp = 1 


(33) 

(34) 

(35) 

(36) 
(37a) 


_9T 

dcp 


dZ 

dW 


[1 + MG(L S (#)] for cp = 1 


(37b) 
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It is convenient at this stage to introduce the new dependent variable © defined by 


© = e x (l-<p) T _ e Ml ~<p) ^ ~ ^ 

A 


(38) 


In terms of this new variable equations (33) to (37b) become 


or 



00 + 3 2 ® - a 2 ® = 0 

dip 2 dcp 2 

(39) 


— = A© at cp = 0 

dcp 

(40) 


1^=0 at ip = 0 

d>p 

(41) 


— = 0 at ip = ip 

dip r 

(42) 

© = 

1 + NH(L (^)) at cp = 1 

(43a) 

^ + A® = 
dcp 

dZ [1 MG(L M))] at cp - 1 

dW L s J 

(43b) 


Equation (39) together with the boundary conditions (40) to (43a) (or alternatively 
(40) to (42) and (43b)) constitute a boundary value problem for © in the rectangular re- 
gion of the W -plane shown in figure 2 which completely determines © (and therefore T) 
as a function of cp and ip. Notice that the particular relation between Z and W (which 
depends only upon the geometry of the porous wall in the physical plane enters only 
through the boundary condition (43a) (or alternatively (43b)) since L g and |dZ/dW| 
must be known functions of ip in order to completely determine these boundary condi- 
tions. It is, however, possible to solve these boundary value problems for arbitrary 
L g and |dZ/dW| (as well as arbitrary H and G), and so the particular geometry and 
boundary conditions can be substituted into the general formulas once this general solu- 
tion has been obtained. 
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General Solution of Boundary Value Problem in W-Plane 


The boundary value problem posed by equation (39) and the boundary conditions (40) 
to (43a) (or alternatively the boundary conditions (40) to (42) and (43b)) can easily be 
solved by the method of separation of variables. To this end we substitute a trial solu- 
tion of the form 


®(ip,cp) =¥(*)*(<?) 


(44) 


into equation (39) and obtain 


ill + = X 2 

^ <f> 

This implies that there exists a constant P such that 


and 


Hence 


$ 


* = c 1 eV^> + c 2 e-' / A^> (45) 

= Cg sin pip + c^ cos Pip (46) 

where Cj to c 4 are arbitrary constants of integration. 

Upon substituting equations (44) to (46) into the boundary condition (40) , we find that 
this boundary condition will be automatically satisfied if 

+ /3 2 - c 2 ^X 2 + /3 2 = X( Cl + c 2 ) 


or 
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II II III III Mill I II mil Mini I II II ■■■II l I in ■ Mil ■ ■■■» 


c 2 " C 1 


j/x 2 + 0 2 -X 
/x 2 + 0 2 + X 


Using this expression in equation (45) gives 
2c* 

$ =- . ± — ; 
f/x 2 + 0 2 + X 

The boundary condition (41) will be automatically satisfied if we put 

c 3 = ° 

and the boundary condition (42) will be satisfied if we set 

0 = — for n = 0, 1, 2. . . . 

^r 

Hence, the solution to both boundary value problems must be of the form 


sinh 


(}fx 2 + 0 2 cp) + |/x 2 + 0 2 cosh ^|/x 2 + 0 2 <p) 


where 



(47) 


(48) 


and the constants C n can be determined so that either the boundary condition (43a) or 
the boundary condition (43b) is satisfied. 


Solution to Boundary Value Problem for Specified Temperature 

First suppose that the boundary condition (43a) applies. Put 
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for n = 0, 1, 2, . . . 


B = CL T X sinh y + y cosh y~\ 
n n|_ 'n 'n 'nl 

Then equation (47) becomes 


© 



X sinh y n cp + y n cosh y n cp 

X sinh v + y cosh v 
'n 'n 'n 



Substituting this result into the boundary condition (43a) yields 


( 49 ) 


1 + NH(L s (i/>)) = 



Hence, the theory of Fourier series shows that 


B 


0 



B 


n 



+ NH(L S ($)] 



for n = 1, 2, 3, . . . 


or 


B Q =1 + NB 0 


B n = NB n for n = 1, 2, 


"N 


> 




where 


(50) 
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B 


1 




h(l s mw 


B 


n 




H(L s (^))cos 



> 


for n = 1, 2, . . . 


( 51 ) 


Substituting equation (50) into equation (49) yields 


© = + N B 

n=0 


X sinh y^cp + y^ cosh y^cp 

X sinh v + v cosh v 
'n 'n 'n 



(52) 


Thus, equation (52) with the constant B n defined in terms of the dimensionless temper- 
ature distribution H by equation (51) is the solution to the boundary value problem. 

In this case, it is of interest for practical applications to have an expression for the 
conduction heat flux q flowing into the surface s . This can be obtained from 

q g = k m n g • vt for (x,y)es 


Hence, for (X, Y)eS (by the use of eqs. (21b), (29), and (30)) 




g-i! =n . vT = V< ? • VT = 


k (t 1 - 1 ) 
m v 1 


|v<p| 


dW 

dZ 


9T 

dtp 


By using equation (38) and the fact that cp = 1 on the boundary s, we find 


w 

dW 

(S® + X0 \ 

k (t, - 1 ) 

m v 1 oo' 

dZ 

W ) 


Upon substituting equations (43a) and (52) into this expression, we obtain 
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q h 
s r 


- 4 J 


dW 


dZ 


<P = 1 


•< 2X + XNH(L (ip))+ N 
s 



n=0 


X cosh v + y sinh y„ 
'n 'n 'n 

X sinh y_ + y„ cosh y„ 
'n 'n 'n 


cosl 


Solution to Boundary Value Problem for Specified Heat Flux 

Now suppose that the boundary condition (43b) applies. Put 

E n = c n + r 2 ) slnh r n + 2Xy n cosh y n] for n = 0, 1 , 2 , . . . 
Then equation (47) becomes 


© 



X sinh y n cp + y n cosh y n <p 
( a2 + y n) 


sinh y n + 2Xy n cosh y R 



Substituting this result into the boundary condition (43b) yields 


dZ 

dW 


[1 + MG(L g 

<P= 1 



Hence, the theory of Fourier series shows that 




Thus, equation (54), with the constants E n defined by equation (55) in terms of the di- 
mensionless heat flux distribution G and the term |dZ/dW (this will be shown by 
eq. (64) to be the reciprocal of the surface velocity), is the solution to this boundary value 
problem. In particular, the upper surface temperature distribution is obtained by using 
equations (54) and (38) evaluated at cp = 1 to obtain 


(*s - { 
Ms 



A. sinh v + y„ cosh v 
'n 'n 'n 

+ ) sinh v + 2Av cosh v 

'n / 'n 'n 'n 



(56) 


HEAT-TRANSFER CHARACTERISTICS FOR A POROUS SQUARE DUCT 

To illustrate the use of the general solution, heat-transfer results will be obtained 
for a square duct with porous walls as shown in figure 3(a) . Cooling fluid is flowing from 
an outer reservoir at temperature t^ through the duct wall and into the duct interior. 
The thickness of the duct wall is used as the reference dimension h r so that a quarter 
of the duct appears in the dimensionless physical plane as shown in figure 3(b) . 


Conformal Mapping Relations Between Physical (Z) and Potential (W) Planes 

The general solution was obtained in the W-plane, and hence, the heat -transfer 
quantities along the duct interior boundary are given in terms of the variable ip. The 
function L g of if) which relates the physical distance in the X,Y-planeto ip must be 
determined. Also the quantity |dw/dz|^_^ must be found as a function of ip since it 
appears in the general solution. These relations can be obtained from the mapping func- 
tion Z — W between the regions in figures 3(b) and 4(a). This mapping has been carried 
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(b) One-quarter of duct in dimensionless physical plane Z = X + iY. 
Figure 3. - Porous cooled square duct. 


(b) Intermediate ^plane u)= + fo?. 



Figure 4. - Mapping of quarter of duct into potential and intermediate planes. 





out in appendix A. The relations between ip and X, Y along the duct interior boundary 
are given by equations (A16). 

The term |dW/dZ | ^ can be found by using the relation 

dW _ dW dco ^ 57 

dZ dw dZ 

It follows from equation (Al) of appendix A that 

WK' - K = 



so that 


™ = L ' . (5, 

dw K ' |/(1 - m 2 )(1~ k 2 cu 2 ) 

The term dw/dZ is given by equation (A2), and the result of combining these relations 
into equation (57) is 


dW _ Vco 
dZ K’Ck 


(59) 


The relation between w and W is given by equation (Al) and the constant C by equa- 
tion (A7). Hence, 


dW 

dZ 



Since W = ip + i on the boundary where w = 1, the identity 


(60) 


sn (ipK' + iK’ - K) = 

ksn (ipK' - K) 

shows that when equation (60) is evaluated on this boundary the following result is ob- 
tained: 
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( 61 ) 


dW = / 2 K ( k 2^ 

dZ o=l \ + k > K ’< k > y isn^K’: 

The absolute value signs are needed under the square root sign because ip < K/K* on 
part of the boundary (between points 6 and 5 in fig. 4(a)). 


Exit Velocity at Inside Surface of Porous Duct 

It follows from equations (13) and (18) that the local velocity in the porous medium 
is given by 


u = iha + j ha 

3X 9Y 

The surface 456 in figure 3(b) is at constant potential, and hence, the exit velocity from 
the porous medium is normal to this surface. It can be seen from the symmetry of the 
problem that the magnitude of the velocity is the same function of distance from the 
corner 5 along the boundary 54 as along 56; so only the former will be computed. Be- 
cause ?(_n/dX - 0 on the boundary 54, it follows that the exit velocity which is in the Y- 
direction is 


v = ha 

S 9Y 


0=1 


(62) 


Since along the boundary 54 


dW 


- i Jilt + i?i£ 


i ha + ^ 

_ dcp 

dZ 

cp- 1 

3Y 9Y_ 

0=1 

ax 3Y_ 

(0=1 


(63) 


it can be seen from equations (62) and (63) that 


V 


s 


dW 

dZ 


l<0= 1 


(64) 


By using equation (6 1) the exit velocity distribution is found to be 
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(65) 


V s = 


{ 


K(k 9 ) 


k(l + k) K’(k) y/sn(i/»K’ - K, k) 


The mapping in equation (A16b) relates ip to X' so that V_ can then be found as a 

* - b S 

function of position along 54. 


Solution With Specified Uniform Temperature on Each Surface of Interior Duct 

Equation (53) provides the dimensionless heat flux along the surface on which the 
wall temperature variation has been specified by equation (7a). To demonstrate the ap- 
plication of this solution, consider the case where one surface of the duct corresponding 
to side 65 j.n figure 3(b) is at the uniform temperature t^ while side 54 is at a different 
uniform temperature tg. Then in the boundary condition (37a) the function H(L g ) of ip 
becomes a unit step function that is equal to zero when ip lies between 0 and K/K'(0 to 
ip / 2) and is equal to unity when ip is between K/K' and 2K/K* ($ r /2 to ip r ). 

The coefficients Bq and B n given by equation (51) now become 

/‘r 

B n = — / d ip = - (66a) 

2 


_ 1 )(n-l)/2 n = 1, 3, 5, . . . (66b) 


Therefore, equation (53) becomes 
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2N 

it 


z 


( _ 1) (n-l)/2 

n 


A cosh v + y sinh y„ 
'n 'n 'n 

A sinh y n + y R cosh y fl 




OO 

A cosh y + y sinh y 
'n f n 'n 

A sinh y + y„ cosh y 
'n n n 

For large n, equation (48) shows that y n ~ mt/ip r and the n n term of the series in 
equation (63) behaves like 


<P 


so that the formal Fourier series in equation (67) does not converge. To remedy this, 
we add and subtract the divergent part of the series to obtain 


r i 

A cos — + — sinh — 

( -l)(n-i)/ 2 ^r ^r CQS ^ 

n ^r A sinh — + IHL cosh — ^r 


(_ 1 } (n-l)/2 JL cos H 
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k m<‘l - * J 


r 


oo 




(-l)( n “l)/ 2 V" 


Vn - — \sinh y n + y n (x - H )cosh y 


n 


X sinh y n + y n cosh y n 


|cos”*f (68) 

K 


The last series in equation (68) now converges, but the series 


OO 



n=l, 3,5,... 


(-D 


(-D/^cos 


(69) 


does not. 

It is shown in appendix B that this series is the formal Fourier series for the func- 
tion (l/2)/cos(77i p/ip ). Hence, when the divergent series is replaced by the singular func- 
tion which it represents, the solution becomes in its final form 


28 




The manipulations carried out above were purely formal since it is not in general proper 
to add and subtract divergent series or to allow them to represent functions. However, 
these manipulations are justified in the case of Fourier series by the use of the theory of 
distributions. This justification is given, for example, in reference 5. 

For the special case when the imposed surface temperature is uniform around the 
entire duct interior, the parameter N in equation (70) is zero. Then equation (70) 
simplifies to 


“s h r 


k (t, - t ) 
m v 1 °° ’ 


dW 


dZ 


2A 


<P= 1 


Therefore, equation (64) shows that for the uniform surface temperature case, 


«s h r 


2Jk (t, - t ) 
irr 1 °° 7 


= V. 


(71) 


Temperature Distribution for Imposed Uniform Heat Flux 

When the imposed heat flux on the inner surface of the duct is uniform, equation (10) 
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shows that M = 0. Consequently, equation (55) becomes in this case 


E 


0 


1 





to W 

dW «=l 


(72a) 


E 


n 






for n = 1, 2, 3, . . . 


(72b) 


In view of the symmetry of the duct, it can be seen that the function |dZ/dW | (which 
is the reciprocal of the exit velocity from the porous material) is an even function of ip 
with respect to the point ip r / 2. The quantity cos(nmp/ip r ), however, is an odd function 
with respect to the point ip r /2 for odd values of n. Hence, E n = 0 for odd values of n. 
Therefore, upon removing the term containing Eq from the sum, the solution for the 
surface temperature distribution given by equation (56) can, in this case, be written as 


h A 


2X 



n=l, 3, 5, 


E_ 


A sum y n - 'n 'n 

(a 2 + y^j sinh y n + 2Xy n cosh y n 


| cos | — ip] (73) 


DISCUSSION 

A general analytical solution was obtained for the heat -transfer behavior in a two- 
dimensional porous medium with coolant being forced through it. The application of the 
solution was demonstrated by examining the heat transfer in a square duct with porous 
cooled walls. Coolant from a reservoir which surrounds the duct is being forced through 
the walls to the duct interior. 

In the general analysis, two of the boundaries of the porous region are at constant 
pressure while the other two have no flow across them. As a result of these conditions 
the porous region occupies a rectangle in the velocity potential plane. The energy equa- 
tion for the temperature distribution in the medium is a separable equation when ex- 
pressed in terms of the independent variables of the potential plane. General solutions 
independent of the specific geometry can then be found in this rectangular region. So- 
lutions of this type were obtained for both the case of an arbitrary specified temperature 
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and the case of a specified heat flux on the porous material boundary through which the 
coolant exits. The general solutions can then be carried from the rectangular region into 
specific physical geometries by conformal mapping. 

Some typical results for the porous cooled square duct which was analyzed in detail 
will now be discussed. First, consider the coolant exit velocity along the duct interior 
surface. This is given by equation (65) as a function of ip, which is related to the phys- 
ical coordinate by equation (A16b). The parameter in equation (65) is k (which also 
determines kj and kg). The constant k is related to the geometric quantity A by 
equation (A8), where A is the ratio of half -width of the outside of the duct to the wall 
thickness. 

The dimensionless exit velocity and duct cross sections for three values of A are 
shown in figure 5. For one -dimensional flow through a plane slab the dimensionless ve- 
locity v uh / /c(p - p ) is unity and the present results approach this limit in the central 

fax U b 

region of the duct side as A becomes large. The velocity is high near the duct interior 
corner because there is less flow resistance in this region as a result of the larger out- 
side surface area and consequently larger cross sectional area for the flow to pass 
through. At the duct interior corner the velocity becomes infinite. As evidenced by the 
deviation of the curve for A = 2. 5 from unity the two-dimensional effects are confined 
to the region within one duct wall thickness of the interior corner. 

Next, consider the heat transfer behavior when there is a specified temperature dis- 



Figure 5. - Dimensionless fluid exit velocity or dimensionless heat flux when surface is at 
uniform temperature. 
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tribution on the inside of the duct. Specifically, it is desired to know for a given tem- 
perature what heat flux occurs at the surface. For a uniform surface temperature, 
equation (71) shows that the dimensionless heat flux q s h r /2Ak m (t^ - t^) is equal to the 
fluid exit velocity. Hence, the ordinate of figure 5 can also be labeled as the dimension- 
less heat flux for the uniform surface temperature case. 

Equation (70) provides results for the situation where adjacent sides of the duct in- 
terior are at different uniform temperatures t^ and tg. For the set of numerical re- 
sults that will be given, the parameter N = (tg - t ^) /(t ^ - t^) was set equal to unity. In 
figure 6, for each value of A there are three curves corresponding to the geometries 
A = 1. 5, 1. 75, and 2. 5. The parameter X = (PCp/2k m ) j«:(p 0 - P s )/mJ regulates the gen- 
eral level of the curves. A large value of p^ - p_ and, hence, large X, is associated 
with a large flow; therefore, the heat flux imposed at the surface is increased for the 
same surface temperature. As in the case of a uniform temperature all around the in- 
terior of the duct, the resulting surface heat flux along each duct side rises toward the 
corner because of the larger flow velocity there. However, very close to the corner there 
is a superimposed local heat conduction effect that alters this trend. Heat flows locally 
by conduction from the side at high temperature t 2 to the side at low temperature t ^ , 
where it must be carried away from the surface in order to maintain the surface temper- 
ature discontinuity. This means that local surface cooling must be provided on the low 
temperature boundary adjacent to the temperature step. This is the region of negative 
heat flux in figure 6. 

Figure 7 gives the surface temperature distribution for the situation where there is 
a uniform heat flux imposed along the entire interior surface of the duct. This solution 
was given by equation (73) and involves the two parameters X and A. The parameter 
A occurs because it determines the parameter k through the relation in equation (A8) 
and k then fixes ip = 2K(k)/K’(k). Results are given for the three values of A shown 
in figure 7 for each of three values of X. 

Since the surface heat flux is uniform, the surface temperature is lowest at the 
corner where according to figure 5 the fluid velocity is high. A larger value of X cor- 
responds to a larger pressure drop and, hence, to greater flow through the porous ma- 
terial. Consequently, the lowest surface temperatures correspond to the largest values 
of X. 
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CONCLUDING REMARKS 


An analytical technique has been developed for determining the heat -transfer behavior 
for a finite two-dimensional porous cooled medium. It was shown that, for the boundary 
conditions imposed herein, every porous material configuration maps into a rectangular 
region in the velocity potential plane. The energy equation is transformed into this plane 
and a general solution obtained. This solution can be mapped into a physical shape by 
using the conformal mapping between the rectangle and the physical plane. 

As a demonstration of the method, heat transfer results were obtained for a square 
duct with porous cooled walls. Coolant is being forced through the walls into the duct 
interior. The two-dimensional effects along the duct interior surface were found to be 
confined to within approximately one wall thickness of the duct interior corner. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, July 17, 1970, 

129-01. 



APPENDIX A 


MAPPING OF QUARTER OF DUCT INTO RECTANGLE OF POTENTIAL PLANE 

In order to relate the solution obtained in the potential plane (W -plane) to the physical 
geometry (Z -plane), the conformal mapping between these two regions must be found. 

The appropriate regions for this geometry are shown in the physical and potential planes 
in figures 3(b) and 4(a). First, it will be convenient to express this mapping in terms of 
the variables in the intermediate co -plane shown in figure 4(b). Here the region occu- 
pies the upper half plane with the boundary along the real (|) axis. Then to aid in some 
of the mathematical manipulation it will be convenient to transform the results into the 
u-plane shown in figure 4(c). In this plane the porous region occupies the entire plane 
with the boundary of the porous region lying along the two cuts on the real axis shown in 
the figure. 

The mapping that transforms the W-plane into the oj - plane in the manner shown by 
figures 4(a) and (b) is given in reference 6 as 

co = sn(WK’ - K) (Al) 

The mapping that transforms the Z -plane into the co-plane in the manner shown by 
figures 3(b) and 4(b) is found by use of the Schwarz -Christoff el transformation to be de- 
termined by 


or 



(A2) 


Z 




dco 

V o>(l - co 2 )(l - k 2 co 2 ) 


(A3) 


where Cj is a constant of integration. 

The integral in equation (A3) is of the type treated by relation 596. 00 on page 268 of 
reference 7. To simplify the integrals, the transformation 
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(o = — (u + %/u 2 - 4k 

2k V 


(A4) 


is used. As shown in reference 7, equation (A3) then transforms into 



(A5) 


The transformation (A4) relates the regions in the o>- and u -planes as shown in fig- 
ures 4(b) and (c). 

The origin of the Z -plane is at point 2, which corresponds to the point «= - ie in 
the u -plane. By integrating from point 2 to any point in the u -plane, equation (A5) 
then provides the following relation between the position in the physical plane and the 
parametric variable u: 



The constant C can be obtained by noting that in figure 3(b) - Zg = i. Equa- 

tion A(6) is then integrated along the branch cut in figure 4(c) from point 3 at u = 1 + k - 
ie to u = 2-y/k and then from u = 2-y/k to u = 1 + k + ie. Equating the result to i gives 
after some simplification 


i = 



da 

^/(a - 2-0c)[(k + l) 2 - a 2 ] 


By use of relation 236. 00 in reference 7, the integral can be converted to an elliptic 
integral. This shows that the constant C is related to the parameter k by 



The quantity k will now be expressed in terms of a physical quantity. The only 
geometric parameter for the square duct configuration is the dimensionless half -width 
A shown in figure 3(b) and equal to Zg - Z 2 - Thus, to obtain an expression for A, 
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equation (A6) is evaluated between points 2 and 3 along the branch cut in figure 4(c). 
This yields 


A = - — 



1+k /*l+k 

da + / da 

|/(« - 2 yfi)[a 2 - (k + 1 ) 2 J |/(a + 2 -\fk)[a 2 _ (k + 1 ) 2 J 


J 


By use of relation 238.00 in reference 7, these integrals are expressed as elliptic in- 
tegrals, and the constant C is eliminated by using equation (A7). When this is done, 
we obtain the following relation between the physical quantity A and the mapping param- 
eter k: 


2 ^) 


_K(k x ) + K(k 2 ) 


(A8) 


where 


k = 1 + k = i - yfc 

1 \/2(l + kj 2 \[2{1 + k) 


In order to complete the solution to the problem, it is necessary to find the relations 

between the values of ip along the boundary cp = 1 and the physical Z -plane. 

For convenience let the coordinates between points 5 and 4 be called X' , Y' while 

s s 

those from 5 to 6 are X g , Y g (fig. 3(b)). Note that Zg is at 1 + i and evaluate equa- 
tion (A6) along the branch cuts in figure 4(c) from point 5 to an arbitrary point between 
either 5 and 4 or 5 and 6. This yields 


Z s 


1 + i + SSJ 

/*oo 

/ d0! 

2 

J a |/(a -2y/k)[a 2 -(k+l) 2 ] 



da 

|/(a + 2yi)[a 2 


- (k + l) 2 ] 


1 + k < a < oo (A9a) 
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Z s = l + i + 



da 


/ 


da 


|/(a - 2y/k)[a 2 - (k + l) 2 ] J a |/(a + 2\/k)[a 2 - (k + l) 2 ] 


-oo<Q!<-(l+k) (A9b) 

By use of 238.00 in reference 7 the integrals are converted to elliptic integrals, the 
constant C is substituted from equation (A7) , and the real and imaginary parts are taken 
to yield 


X s = 1 


Y s= 1 + 


- — F /sin -1 k \ 

( k 2 ) y y -a + (1 + k) y 

in 

y -a + (1 + 


S. -oo < a < -(1 + k) (A10a) 


k) 


X s = 1 + 


2K(k 2 ) 





-F 



2(1 + k) 
a + (1 + k)’ 



> 1 + k < a < 


oo 


(AlOb) 


y ; =i J 

Now it follows from equation (A4) that 2ku> - u = |/u 2 - 4k. Square both sides and 
solve for u to obtain 


u = ku> + — 


Ci> 


(All) 
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In the ranges 1 + k < a < «> and -°° < a < -(1 + k), for which equation (A10) applies, 
the points are on the real axis of both the u- and oj - planes so that equation (All) gives 

a = k£ + A (A12) 

This is substituted into equation (A 10) to give 




Equation (Al) relates the oo- and W-planes and can be used to make the final 

transformation of equation (A13) and relate X , Y and X' , Y' to ip. For X . Y_ 

s s s s s s 

the coordinates are between points 6 and 5, and from figure 4(a), W = ip + i where 0 ^ 
ip < K/K’. Then 

| = sn(ipK' + iK T - K) 
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which gives 




1 1 
k sn(K - tpK') 


0 < ip < 


-OO < £ 


JC 

K’ 


1 

k 


Similarly, between points 5 and 4 




k sn(^K' - K) 


' I^<2l 

K' K’ 


— < 2j < oo 

k 


Equations (A14) are used in (A13) to eliminate Note that 


-.1/2 r* 


2(1 + k) 


±k£ ± — + 1 + k 


-.1/2 


2_ 

±£ 


1 + k 




Substitute either (A14a) or (A14b) and define S2 as 


(A14a) 


(A14b) 


fi s. 


2(1 + k/lsn^K’ - K,k) 


1/2 


£l + k | sn(i/)K' - K, k) | J 1 + | sn^K’ - K, k) | J 


Then equations (A13) become 


X s = 1 


Y s = 1 + 


2K(k 2 ) 



J 


(A15) 


(A16a) 
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(A16b) 



42 






If we cross multiply, 

(z - e -aT )y(z) = k(z - l)x(z) 

divide by z, 

(* ■ ° k ( x ' l) x(z) 

transform to the discrete time domain, remembering that (1 /z)y n = y n _i 

_ oT 

y n - e y n -i = k ^ x n - x n-i) 

and rearranging, the final form of the difference equation is obtained: 

Yn = e" aT y n -i + k Cx n - x n _j) 

At the sampling points, the difference equation is the exact solution for the response of the 
equivalent analog system to a stair-step function input. Therefore, the solution does not become 
numerically unstable for large sampling intervals. However, the output of the system will be realistic 
only if the sampling is done frequently enough to result in a stair-step waveform which is a good 
approximation of the original continuous input waveform. 

For comparison, the equations for the above filter using numerical integration (fourth-order 
Adams-Bashford) are: 

yn = yn-l + 24 - ^9yn-i + 37y n _ 2 - 9y n -3) 


y n = -ay n + kx n 

The computation efficiency of the difference equation is obviously better than that of the 
numerical integration technique. Also, this numerical integration method becomes numerically 
unstable for large sampling intervals, as confirmed during the investigation. 

Difference equations were used to simulate the aircraft control and automatic stabilization 
systems, to replace the analog autopilot with a digital version, and to generate the turbulent wind 
components. 

The aircraft control system simulator consists of a throttle servo with engine lag and an 
elevator servo with an elevator surface lag. Both of these servos have the same general form, a 
second-order transfer function with 0.7 critical damping in series with a first-order low pass filter. 
The form of the total transfer function and the corresponding difference equation are shown in 
section 7 of table 8. 

The washout filter in the stability augmentation system discussed above is simulated by the 
difference equation of section 1 in table 8. 
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_ J. 

1 + exp /- jjL + iW - icj -—j. 1 + exp <f- JjL + iW + icj — .'j, 


= 2iJ(-l)-e' (n,r/ ^ )(1 “ <P) e -^r)^ sin n,c 
n^O 




for cp < 1 (B2) 


Upon combining the two fractions on the left side of equation (B2) and collecting 
terms, we find 


1 


+ exp,|- JjL + iW - icj — ^ 1 + exp <f- jjl + iW + icj — 


= - 1 sin 


7 TC 


fcos 

— (ip + icp - i) 

+ cos 

l 

A J 



Hence, 
_ 1 

equation (B2) becomes 

sin — 

^r 

2 cos 

— (ip + icp - i)l + cos — 


A J ^r 


OO 

= £<- 


n=0 

Or upon taking real parts, 


-(mr/ip r )(l-(p) -(niriip)/ip 1 


e - e r sinl nnc 


for cp < 1 
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2 


Re 


sin 


7TC 


7TC 

COS — + COS 


^r 


J~(ip + i(p - i) 

L*r 


QO 

I 

n=0 


(-D n 


-(nn/ip )(l-(p) 


sin !HL£cosi^ for <p < 1 




Upon taking the limit cp — 1, we find that relation (Bl) holds. Of course this limit must 
be taken in the sense of distributions, in which case the interchange of the limit assoc- 
iated with the summation and the limit associated with cp is justified. 
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complete, and a lasting contribution to existing 
knowledge. 

TECHNICAL NOTES: Information less broad 
in scope but nevertheless of importance as a 
contribution to existing knowledge. 

TECHNICAL MEMORANDUMS: 
Information receiving limited distribution 
because of preliminary data, security classifica- 
tion, or other reasons. 

CONTRACTOR REPORTS: Scientific and 
technical information generated under a NASA 
contract or grant and considered an important 
contribution to existing knowledge. 


TECHNICAL TRANSLATIONS: Information 
published in a foreign language considered 
to merit NASA distribution in English. 

SPECIAL PUBLICATIONS: Information 
derived from or of value to NASA activities. 
Publications include conference proceedings, 
monographs, data compilations, handbooks, 
sourcebooks, and special bibliographies. 

TECHNOLOGY UTILIZATION 
PUBLICATIONS: Information on technology 
used by NASA that may be of particular 
interest in commercial and other non-aerospace 
applications. Publications include Tech Briefs, 
Technology Utilization Reports and Notes, 
and Technology Surveys. 


Details on the availability of these publications may be obtained from: 

SCIENTIFIC AND TECHNICAL INFORMATION DIVISION 

NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 

Washington, D.C. 20546 


